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ABSTRACT 


Ttiis report describes an algorithm to carry out sensitivities, 
uncertainties and overall imprecision studies to a set of input 
parameters for a one dimensional steady state ozone 0iotocheinistry 
model. 'Ibis algorithm can be used to evaluate steady state 
perturbations due to point source or distributed ejection of H20, CLX, 
and NOX, besides, varying the incident solar flux. liiis algorithm is 
operational on IBM OS/360-91 computer at NASA/Goddard Space Flight 
Center's Science and Applications Computer Center (SACC) . 
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1. INTRODUCTION 


Rundel, Butler and Stolarski (1) have carried out stratosj^eric 
.ozone uncertainties studies using a Monte-Carlo method. Hiis algorithm 
uses an atmosfiiere covering altitudes 0-60 km for user specified 
thicknesses of 1.25 to 5.0 km. liie algorithm incorporates odd oxygen, 
nitrogen, hydrogen and chlorine chemistry and utilitizes 58 chemical 
reactions described in Table 1. Ttie fixed (user input) profiles of N2, 

02 and C02 (stable molecules) are used. The concentrations of N20, H20, 

H2 and CH4 are determined from self consistent diffusion calculations 
with fixed boundary conditions for the above species. The OX, 
and CO are calculated by solving differential equations using boundary 
conditions chosen to reproduce the atmospiieric measurements. 

For analyzing the uncertainties of stratospheric perturbations, 
this algorithm uses techniques of decoupling iterations. The 
differential equations solved are; 

a. diffusion equations 

b. local f^ioto-chemical equations for species or ratio of 
species. 

This algorithm is capable of carrying out uncertainties studies for 
an altitude range of 0-60 km in steps of 1.25 km. Table 1 shows the 
basic chemistry of the odd oxygen, nitrogen, hydrogen, chlorine compounds 
together with the methane oxidation chain. 

In general the steady state distribution of an atmospheric trace gas 
(j) is calculated by solving the continuity equation: 

V 4>j, == 

where is he number density of the jth specie, the flux vector, 

Pj and Lj are the chemical production and loss terms. SOi^enotes the set of 
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TABLE 1 


REACTIONS USED IN OUSA ALGORITHM 



Reactions 

Reaction Number 

1 . 

Cl + 0^ ^ CIO + 0 ^ 

(42) 

2 . 

CIO + 0 ^ Cl + O 2 

(43) 

3. 

CIO + NO -»■ Cl + NO 2 

(46) 

4. 

Cl + H 2 HCL + H 

(40) 

5. 

Cl + N02^ HCl + O 2 

(41) 

6 . 

Cl + CH^ HCl + CH^ 

(44) 

7. 

Cl + H 2 O 2 ^ HCl + HO 2 

(13) 

8 . 

HCl + OH ^ Cl + H 2 O 

(45) 

9. 

CIO + NO 2 + CIONO 2 + M 

(52) 

10 . 

CIONO 2 +0^ CIO + NO^ 


11 . 

CIONO 2 + OH ^ HOCl + NO^ 

(25) 

12 . 

0 + O 2 + M 0^ + M 

( 2 ) 

13. 

0 + 0 ^ O 2 + O 2 

( 4) 

14. 

NO + 0^ ^ NO 2 + ^2 

(15) 

15. 

NO 2 t O 3 - NO 3 + O 2 

(16) 

16. 

OH + O 3 ^ HO 2 + O 2 

( 9) 

17. 

HO 2 + 03 ^ OH + O 2 + O 2 

( 10 ) 

18. 

H + O 3 HO + O 2 

( 22 ) 

19. 

NO 2 +0 NO + O 2 

(17) 

20 . 

NO 2 + 0 + M ^ NO 3 + M 

(19) 

21 . 

OH + 0 H + O 2 

( 7) 

22 . 

HO 2 + 0 OH + O 2 

( 11 ) 

23. 


(55) 

24. 

HNO 3 + 0 OH + NO 3 

(48) 

25. 

CH 2 O + 0 HCO + OH 

(37) 

26. 

27. 

28. 

N 2 O + 0(^D) ^ NO + NO 
H 2 O + 0 (^D) OH + OH 

H2 + 0(^D) -»-0H + OH 

( 1 ) 
( 5) 
(26) 
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29. 

30. 

31. 

CH^ + 0(^D) ^ OH + H 
CH^ + 0(^D) -y + CH^O 

0(^D) + M 0(^P) + M 

(21) 

(54) 

(14) 

32. 

NO + NO^ ^ NO 2 + NO 2 

(49) 

33. 

N + NO N 2 + 0 

(35) 

34 . 

N + O 2 ^ NO + 0 

(36) 

35. 

N + 0^ NO + O 2 

(56) 

36. 

NO + HO 2 -> NO 2 + OH 

( 3) 

37. 

NO^ + OH + M ^ HNO^ + M 

(18) 

38. 

OH + HNO^ -S' NO^ + H 2 O 

(20) 

39. 

OH + H 2 O 2 ^ HO 2 + H 2 O : 

(47) 

40. 

OH + HO 2 ^ O 2 + H 2 O 

( 6) 

41 . 

OH + OH 0 + H 2 O 

(12) 

42. 

HO + H 2 -V H + H 2 O 

(50) 

43. 

OH + CH, CH- + H»0 

4 3 2 

(29) 

44. 

OH + CH 2 O ^ HCO + H 2 O 

(38) 

45. 

OH + CH^Cl -► Ch2Cl +^.H20 

(32) 

46 . 

OH + CO ->-^-C02 + H 

(28) 

47. 

OH + OH + M -> H 2 O 2 + M 

(23) 

48. 

H + O 2 + M HO 2 + M 

(8) 

49. 

HO 2 + H02-»- H 2 O 2 + O 2 

(24) 

50. 

CH^ + 0^ + M -*■ CH^O^ + M 

(30) 

51. 

CH^O^ + NO ^ CH^O + NO 2 

(31) 

52. 

CH 2 O 2 + NO 2 CH^O + NO^ 

(53) 

53. 

CH 2 O 2 + HO 2 -»■ CH^OOH + NO 2 

(33) 

54 . 

CH^O + O 2 CH 2 O + HO 2 

(34) 

55. 

HCO + O 2 ->CO + HO 2 

(39) 

56. 

CH^ + 0(^D) ^ CH^ + OH 

(27) 

57. 

NO 2 + NO^ + M ^ N20^ + M 

(21) 

58. 

N20^ + M -> NO 2 + NOj +"M 

(51 
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all reactants v^ich enter into the chemistry of the jth specie. The net 
upward eddy flux is approximated by , . 

v^ere is the vertical eddy coefficient. [M] is the density of air and-fj,= 
denotes the volume mixing ratio of a given constituent , 

2. C0MPUTATI®4 PROCEDURE 

Figure 1 illustrates the flow-chart of the calculations of ambient 
constituents concentrations in the model and the method used for chemical 
families and obtaining self-consistent finally conveged species by 
iteration. The calculations sequence consists of: 

1. Initial (starting) values of 03, BOX, NOX, CLX etc. 

2. Calculate column content of 03 

3. Calculate J-coefficients 

4. Calculate 0/03, 0(lD)/03, 03/0X (f^oto-chemical equilibrium) 

5. Calculate new OX (diffusion) 

6. Calculate N20 (diffusion) 

7. Calculate BOX (OB, B, B02, B2, 02 photo-chemical equilbrium) 

8. Calculate CB4, B2, B20 (diffusion) 

9. Calculate CB4 - oxidation products (photo-chemical 
equilbrium) 

10. Calculate CO (diffusion) 

11. Calculate CC14, CB3C1, CF2C12, CFCL3 (diffusion) 

12. Calculate CLX (diffusion) 

13. Calculate CL, CLO, BCL, CLCW02, from CLX (photo-chanical 
equilibrium) 

14. Calculate odd nitrogen (diffusion) 

15. Calculate NO, N02, N03, BN03 from odd nitrogen 
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16. Recycle to 2. 

Imprecision of rate constants used for the reactions in Table 1 may 
range from small (10 to 15%) to large (up to factor of 10) variation from 
the mean values. Such a large degree of variation is most coveniently 
dealt with by using logarithmic variation. Monte-Carlo calculation 
necessitates a probability function for each reaction rate, OUSA 
algorithm uses a Gaussian distribution in log-k (reaction rate) and the 
probability distribution is given by „ 

( 3 ) 

where e- is the standard deviation. 

2.2 Subroutine Main (OUSA) 

This subroutine is a driver program calling most of the subroutine 
detailed in Figure 2, The subroutines called by this routine are: 


1. 

Images 

12. 

HOX 

2. 

ERRSET* 

13. 

CH4D 

3. 

INPUT 

14. 

CH4 

4. 

DIFFUS 

15. 

CLX* 

5. 

OUTl 

16. 

NOX 

6. 

GGNRF 

17. 

CVAL 

7. 

COMPK 

18. 

OUT2 

8. 

03C0L 

19. 

OUTPUT 

9. 

JCOEFF 

20. 

RDPERT (FLUXBC, POINT 




DIST, SOLAR, MIXING) 

10. 

OX 

• 

CM 

OUT3 

11. 

N20D 




♦ERRSET - is system: routine to prevent printout due to underflows. 
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Figure 2 - Block Diagram For OUSA Algorithm 
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Figure -2 Block Diagram For OUSA Algorithm 
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Besides calling for input and output operations via subroutines; 
INPUT, OUTl, (XJT2, 0UT3, and CXJTPUT, this program initiates the routines 
to provide normal distribution of random number (DO L(X)P 111) . It 
provides the initial values of various species, productions parameters, 
reactions dissociation rates etc,,. 

Tbe statements 16300-18100 carryout computations for the ambient 
atmosphere and the convergence of ozone amounts at each atmospAieric 
levels. Depending upon the input print out criterion, the various 
species concentration and mixing ratios are printed and the uncertainity 
studies species are stored on a magnetic tape. 

Ttie statements 20700-31200 repeat the computations identical to the 
ambient atmosphere (detailed in the preceding paragraph) for a perturbed 
atmosphere by using the following selection of perturbations, invoked via 
the subroutine RDPERT; 

1. Flux boundary 

2. Point Sources 

3. Distributed Sources 

4. Variation of Solar Irradiance 

5. Variations of Mixing Ratios. 

Hie unformatted outpnjt of this program is stored on the magnetic tape 
(unit-8) and is used with the program STATIS for canpxiting the statistics 
of Monte-Carlo simulations of ozone uncertaintity studies. 

2.3 Subroutine IMAGES 

This subroutine is onployed to obtain the echo-check of the input 
data specified for a given data unit. For the OUSA algorithm the input 
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resides on unit 5 and 11. It consists of initial atmosgiieric parameters 
and perturbations related parameters. 

2.4 Subroutine INPUT 

Ibis subroutine provides the initial atmospheric parameters for the 
execution of OUSA algorithm. Table 2 shows the INPUT parameters together 
with their explanations. 

2.5 Subroutine DIFFUS 

The main purpose of this subroutine is to ccnipute the diffusion 
related parameters and the molecular number density at all the levels of 
the atmosphere using Mo (number density at bottom) values. 

2.6 Subroutine COMPK 

This subroutine is used for calculating two body and three body 
reaction rates in accordance with a NASA report (2) . Depending on the 
algorithm central line option or uncertainty (Monte-Carlo) simulations, 
these reaction rates are multiplied by a log-normal random distribution. 

2.7 Subroutine OUTl 

This subroutine is employed to print the initial atmospheric 
constituents, as a function of altitude and parameters related to "run" 
setup, the reaction rates related variables, diffusion parameters etc.,. 

2.8 Subroutine 03C0L 

This subroutine is employed for calculating the columnar ozone 
amounts at tlie various levels of the atmosphere. 

2.9 Subroutine JCOEFF 

Ihe photo-dissociation coefficients for various reactions shown in 
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TABLE 2 


INPUT PARAMETERS 


1 . NOPERT 

2. IRAND 

3i IMIX 

4. IFLX 

5. 2BOT 

6. ZINC 

7 . NZ 

8. M(l) 

9 . NTROP 

10. TT(1) 

11. Kl, Al, Bl, K3, K4, K2 

12. VOLN 

13. FLXFAC 

14. TEST 

15. FBC(I) 


-Number of kind of perturbations 
-Number of Monte Carlo simulation 
for each perturbation 
-Flag for printing of mixing ratio 
boundary i.e., if IMIX=0 print 
-Flag for printing for flux boundary 
conditions i.e., if IFLX=0, then 
print 

-Altitude of bottom of atmospiiere 
-Altitude increment (in km) for each 
level of atmosphere 
-Number of atmospheric levels 
-Molecular number density at the 
bottom of the atmosphere 
-Level number for the troposheric 
layer 

-Input temperature at various 
levels 

-These are eddy diffusion 
parameters 

-Diffusion related parameter 
-Flux factor 

-Convergence test value for the 
ozone amount at each level 
-Boundary values of fluxes 
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MR(I) 


-Boundary values of mixing ratios 

17. lOPT(I) -Options for flux and mixing ratio 

boundary conditions 

18. QCOI(I), QCH4(I) ,QH2(I) -Values of carbon monoxide, 
QH20(I), QCLX(I), QN20(I) methane, hydrogen, water, sun 

chlorine specie, nitricoxide, 

19. RAIN(I), SOL(I) -Rain out and solubility para- 

meters as a function of height 

20. SI (I) . . . .S17 (I) -Absorption cross sections and 

J(1,I) . . . .J(22,I) rates as a funtion of height 

21. K(I), E(I), 1=1,60 -Rates and energy parameters for 

reactions 

22. U(I), USEN(I), USENH(I); -Reactions rates uncertainty 

1=1,60 parameters 
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Table 1 are calculated using the solar flux and absorptions cross- 

O 

section. ‘Ibe solar spectrum is divided into 25 regions beyond 3200A 
These calculations are done by assuming no attenutions of incident solar 
radiation. The incident solar flux below 3200 A is assumed to be 
absorbed by 02 and 03 and the attenuation due to other species is assumed to be 
optically thin. The {*ioto-disociation coefficients are diurnally 
averaged by using Rundel's technique (3). 

2.10 Subroutine OX 

This subroutine computes the daytime and 24 hour averages for 

concentration of 0, 0(1D) , and 03. After each entry to this subroutine, 
the loss and production terms (CCi and DDi) are computed and upper and 

lower boundary conditions are set. The diffusion computations are then 

accomplished by using subroutine DIF. 

2.11 Subroutine DIFCAL 

This subroutine is used for computing the steady state values of 
various species by taking into account eddy diffusion. This routine has 
following entry points to carry out diffusion comp>utation of various 
species: 



ENTRY POINTS 

PURPOSE 

1. 

COD 

Carbon monoxide 

2. 

N20D 

Nitrous oxide 

3. 

NOXD 

Odd nitrogen 

4. 

CH4D 

Methane 

5. 

CLXSOR 

Carbon tetra-chloride 

6. 

CLXPER 

CFC13 and CF2C12 
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7, CLXDIF Odd chlorine 

8. OXD Odd oxygen 

The calculations for diffusion of hydrogen and water are carried out 
after the completion of methane diffusion. 

After each entry to the above points the loss and production terms 
(CCi and DDi) are computed and upper and lower boundary conditions are 
set. The diffusion computations are then accomplished by using 
subroutine DIF. 

2.12 Subroutine DIF 

Ihis subroutine solves the differential equation by using the 

Nicholson method. Depending on the upper and lower boandary conditions for 

each specie, the differential equation is then solved for the all the 
atmosfieric levels. Tiie variable Ri (i=l, NZ) returns the mixing ratio 

values for the specie desired in diffusion computation. 

2.13 Subroutine CH4 

This subroutine is used for obtaining the values of CH3, CH302, 
H2C0, and HCO species. In order to carry out diffusion computation for 
the methyl hydrogen peroxide (CH3COOH) and for carbon monoxide, this 
subroutine makes call to entries MOOHD and COD via subroutine DlfWIS. 

2.14 Subroutine NOX 

This subroutine carries out the NOX (N, NO, N02, N03, N205, CL0N02 
and HN03) computation using local photo-equilibrium and CLCW02/N02, 
HN03/N02, N03/N02, N0/N02, and N02/odd nitrogen ratio computations. 

N205 day-day and night— night effects are also used in these computations. 
This subroutine also calls for the entry to DIF program to perform 
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diffusion calculations for the odd nitrogen. 

2.15 Subroutine HOX 

This subroutine carries out the HOX (OH + H02 + HCO + H) related 
computations using steady state values and day-night effects. Hiis 
subroutine also invokes the entry point H202D to DIFMIS subroutine where 
the hydrogen peroxide diffusion calculations are carried out. 

2.16 Subroutine CLX 

This subroutine carries out source calculations for CLX solving 
differential equations for each source molecule assuming a fixed mixing 
ratio at specified height and zero flux at the top of the atmosf^ere. 
This subroutine calls for the entry points CLXSCffl (selection for CC14 and 
CH3C1) and CLXPER (for Flouro-Carbon-11 and Flouro-Carbon 12) diffusion 
calculations in subroutine DIFCAL. This subroutine accounts for the day 
effects of HOCL/CL, CLO/CL and day night effects for CLCN02 and CLO. It 
is assumed that HCL concentrations do not vary on a 24 hour time scale. 
The day is approximated by 12 hours of constant insolation and night as 
12 hours of complete darkness. It is also assumed that CL/CLO approaches 
a rapid equilibrium. 

2.17 Subroutine DIF2 

Tiiis subroutine is used for calculating diffusion of water. 

2.18 Subroutine DIFMIS 

This subroutine is used for calculating the hydrogen 
peroxide and methyl hydrogen peroxide diffusion. 

2.19 Subroutine 0UT2 

This subroutine prints out the number density and mixing ratio as' a 
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fmction of height for the converged ambient atmosf^ere iterations. 

2.20 Subroutine OUTPUT 

This subroutine carries out OX loss and production terms and prints 
them as a function of height. 

2.21 Subroutine CVAL 

Ihis subroutine carries out mixing ratio calculation for a large 
number of molecular species. 

2.22 Subroutine RAINCXJ 

This subroutine computes rainout concentrations of water, odd N, 
H202 and odd CL. 

2.23 Subroutine 0UT3 

This subroutine prints out the densities and mixing ratios for various 
constituents for converged perturbed atmosj^ere. 

2.24 Subroutine RDPERT 

Hiis subroutine provides user-option perturbation parameters input 
for flux boundary conditions, injection of various specie at different 
altitudes, specie profile changes, solar insolation variation, and mixing 
ratio boundary parameters alterations. 

2.25 Subroutine PERTUR 

Ihis subroutine enters the user specified perturbations input 
parameters. Ihe values of these input parameters are printed out as well. 

2.26 Subroutine GGNRX 

Ihis is an IBM OS/360 assembly language program to obtain fast 
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normal random number generation. It calls for assembler control section 
GGNRF (N,A) , where A sets the seed for N normal random number generation. 
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4. PROGRAM JOB COJTROL LANGUAGE (JCL) CARDS 

The following listing details the IBM OS/360 Job Control Language 
(JCL) cards needed to successfully compete the OUSA program. Various 
subroutines detailed in section 2 are preserved under the library ZIOPB 
LIB 0M03. Ihe input data is stored on member DATA3 and perturbation 
parameters are accessible through the member PERTZERO. The outpot of 
this run are stored on a magnetic tape for future statistical analysis. 


A/ZIOPBRUN JOB ( L A004 1 1 95 A , T » OPBOO 1 t HO 1 H03 > t 059 , MSGL E VEL =1 00 
//*C0MP-RUN ' 00 
/'/FORT EXEC F0RTRANH,PARM='MAP»XREF* 00' 
//SOURCE. SYSIN DO D I SP = SHR .OSN = Z lOPB ,L I B .0M0^3 ( DM A I N > - 00 
H DD DISP=SHR,0SN=Z10PB.LIB.0M03(DIMAGES) 00 
// OD D1SP=SHR.0SN=Z10PB.LIB.0M03(DIN) ^ qq; 
H DD 0ISP = SHR,DSN = Z10PB.LIB.0M03(DDIFFUS ) 00^ 
// DD DISP=SHR,DSN = Z10PB.LIB.0M03(DC0MPK ) 00! 
H DD DISP = SHR,DSN = Z10P8.LIB.OMO3(DOUTl ) 00< 
// DD DISP=SHR,DSN=Z10PB.LIB.0M03(D03C0L) 00' 
H DD DISP=SHR,0SN=Z10PB.LIB.0M03(D0UAD) 00' 
// OD DISP=SHR.DSN=Z10PB.LIB.OM03(DJCOEFF) OOf 
H DD DISP=SHR,DSN=Z10PB.L IB.0M03(D0X) 00' 
// DD DISP=SHR,DSN=Z10PB.LIB.0M03(DDIFCAL) OOf 
/'/ DD DI SP = SHR.DSN = Z10P8.L IB.0M03(DDIF ) 00- 
// DD 0ISP=SHR,0SN=Z10PB.LIB.0M03(DCH4) 00' 
H OD DISP=SHR,DSN=Z10PR.LIB.0M03(DN0X) . 00' 
// OD 0ISP=SHR,0SN = Z10PB.LI6.0M03(DH0X ) 00 
H DD DISP=SHR,0SN=Z10PB.LIB.0M03(DCLX) 00' 
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// DD DISP=SHR«DSN=Z10PR.LIB.0M03(0niF2) 00 
H DD DISP=SHR,DSN=Z10PB.L IB.0M03(DDIFMIS ) 00 
// DD DISP=SHR.DSN=Z10PB,LIB.OMO3(D0UT2) 00 
H DD DISP = SHR,DSN = Z10PB.L IB.0M03(D0UTPIJT ) 00 
// DD DISP=SHR,DSN=Z10PB.L1B.0M03(DCVAL) • 00 
H DD DISP=SHR,DSN=Z10PB.LIB.0M03(DRAIN) . 00 
// DO DISP=SHR,DSN=Z10PB.LIB.0M03<D0UT3) 00 
H DD DISP=SHR,DSN=210PB.LIB.0M03(ORDPERT) 00 
// DD DISP=SHRtDSN=Z10PB.LIB.0M03(DPERTUR) 00 
/♦ 00 
//ASMG EXEC ASMG 00 
/•/SOURCE. SYSIN DD D I SP = SHR , DSN= Z lOPB . L I B .0M03 ( GGNRF > 00 
/♦ , 00 
/7TR0P0 EXEC LOADER , PROG= L I NKER . P ARM= • S I ZE = 380K ♦ E P = MA I N • t R EG I ON . G0 = 380K 00 
//GO.SYSLIB DO D I S P= SHR , DSN=SYS2 . I MSLS 00 
/^/GO.FTOSFOOl DD 0 1 SP = SHR , DSN = Z 10P8 .L I B . 0M03 ( DATA 3 ) . L ABEL = ( » , , I N ) 00 
//■GO.SYSUDUMP DD SYSOUT = A .SPACE = (CYL ♦( 2,2 ) » ' 00 
//G0.FT08F001 DO DSN= 2 1 OPB. 8 AHE TH I . D I SP= ( NE W , KE EP ) . UN I T= ( 6250 , . DE F ER ) , 00 
// LABEL=(A,NL,,OUT) . DCB = ( REC FM=VBS ♦ LRECL=32AA . BLK S I Z E= 1 622A , 0 EN= 3 ) ♦ 00 
H V0L=SER=DB0002 00 
//GO.FTl IFOOl DD D I S P= SHR , DSN = Z lOPB . L I B . D AT A ( PER T Z ERO ) ♦ L ABE L = ( , , * I N ) 00 
/-« 00 
//NOTIFY EXEC NOT I F YT S « MSG= • TROP DONE* 00 
H 00 


5. FORTRAN LISTING OF ALGORITHM 

SEE THE ENCLOSED COMPUTER PRINTOUT 

6. INPUT PARAMETERS 

SEE THE ENCLOSED COMPUTER PRINTOUT 

7. OUTPUT FROM TEST RUN 


SEE THE ENCLOSED COMPUTER PRINTOUT 


